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Spatial and Modal Superconvergence of the Discontinuous 
Galerkin Method for Linear Equations 

N. Chalmers*and L. Krivodonova^ 


Abstract 

We apply the discontinuous Galerkin finite element method with a degree p polynomial 
basis to the linear advection equation and derive a PDE which the numerical solution solves 
exactly. We use a Fourier approach to derive polynomial solutions to this PDE and show that 
the polynomials are closely related to the Pade approximant of the exponential function. 
We show that for a uniform mesh of N elements there exist (p + 1)N independent polynomial 
solutions, N of which can be viewed as physical and pN as non-physical. We show that the 
accumulation error of the physical mode is of order 2p + 1. In contrast, the non-physical 
modes are damped out exponentially quickly. We use these results to present a simple proof 
of the superconvergence of the DG method on uniform grids as well as show a connection 
between spatial super convergence and the superaccuracies in dissipation and dispersion errors 
of the scheme. Finally, we show that for a class of initial projections on a uniform mesh, the 
superconvergent points of the numerical error tend exponentially quickly towards the downwind 
based Radau points. 


1 Introduction 

In this paper we apply the discontinuous Galerkin (DG) method with a polynomial basis to the 
one-dimensional linear hyperbolic problem 


ut + aUx = Q ( 1 ) 

with a > 0 constant, subject to periodic boundary conditions on interval / and sufficiently smooth 
initial data rt(a:,0). We derive a partial differential equation (PDE) on the j-th cell which is solved 
by the polynomial numerical solution Uj exactly. This PDE is similar to the original advection 
equation with a forcing term. Then, by applying classical Fourier analysis, we hud polynomial 
solutions to this PDE. These solutions are closely related to both the {p + l)-th right Radau 
polynomial Rp+i and the Pade approximant of the exponential function e^, where p is the 
degree of the polynomial approximation. We use these particular polynomial solutions to show 
that for a uniform computational mesh of N elements there exist {p+ 1)N independent polynomial 
solutions, N of which can be seen as physical and pN as non-physical. Furthermore, we prove a 
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property of the super convergent points of the DG method conjectured by Biswas et al [5] in 1994. 
Specihcally, we show that for a class of initial projections on a uniform mesh the numerical error 
will tend to a superconvergent form with order p + 2 at the right Radau points. Moreover, we 
establish conditions on the projection of the initial data for which this property holds. We also 
show that this superconvergent form has high-order accuracy in the moments of the solution, i.e., 
the projection of the numerical solution onto the m-th Legendre polynomial is order 2p -|- 1 — m 
accurate. In particular, the cell average of the solution is advected with order 2p -|- 1 accuracy. 

Superconvergence of the DG method for one-dimensional problems has been studied in several 
papers. Following the conjecture made by Biswas et al, Adjerid et al in [2] proved order (p -|- 2) 
convergence of the DG solution at the downwind-based Radau points and order 2p-|- 1 convergence 
at the downwind end of each cell for ODEs. An order (p + |) convergence rate of the DG solution 
to a particular projection of the exact solution was later shown by Gheng and Shu in [HI [7]. Yang 
and Shu then showed the same super convergence property with order p -|- 2 convergence for linear 
hyperbolic equations in [13]. Fourier analysis of the DG solution has also been applied to the DG 
solution in order to investigate superconvergence by symbolically manipulating the discretetization 
matrices for low order (p = 1,2, and 3) approximations [HI |8]. 

Likewise, the connection between the DG scheme and the Fade approximants of the exponential 
function has been observed in several works. The stability region of the DG method for ODEs 
was demonstrated by Le Saint and Raviart [12] to be given by |R(A/i)| ^ 1 where R{z) is the 
Fade approximant of and h is the grid spacing. In [9], Hu and Atkins conjectured that certain 
polynomials involved in the analysis of the numerical dispersion relation are related to Fade 
approximant of exp( 2 ;) and used this to show that the numerical dispersion relation is accurate 
to where k is the wavenumber. This conjecture was proven and an extended analysis of 

the dispersion and dissipation errors was given by Ainsworth in [3]. Later, a connection between 
the spectrum of the DG method on linear problems and the Fade approximant of the exponential 
was investigated by Krivodonova in Qin in [10]. In this paper, we demonstrate that the numerical 
solutions of the DG method are themselves closely related to this Fade approximant and furthermore 
both the superconvergent local errors and superaccurate errors in dissipation and dispersion of the 
method can be seen as resulting from the accuracy of this Fade approximant. We also use these 
numerical solutions to present a new proof of spatial super convergence of the DG method on uniform 
grids. This proof is simpler than existing ones and does not require the assumption that the 
initial projection is superconvergent at the beginning of the computation. Moreover, we show what 
condition the initial projection should satisfy for the error to become and remain superconvergent 
at the right Radau points after sufficient time. 

The remainder of this paper is organized as follows. In Section 2 we begin by deriving the DG 
method and use the formulation to write a FDE which is solved by the numerical solution Uj on 
the j-t\i cell. The main body of this paper is in Section 3 where by decomposing the numerical 
solution into Fourier modes, we use this FDE to hnd particular numerical solutions and establish 
our main results. These results are then illustrated numerically in Section 4. Finally, conclusions 
and thoughts on future work are discussed in Section 5. 
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2 DG Method 


Our goal in this section is to derive a partial differential equation for the numerical solution of the 
DG method. We begin with the linear conservation law ([T]). The domain is discretized into mesh 
elements Ij = [xj,Xj^]\ of size hj = Xj+i — Xj, j = 1,2, The discontinuous Galerkin spatial 

discretization on cell Ij is obtained by approximating the exact solution u by a polynomial of degree 
p, Uj. Multiplying ([T]) by a test function V EVp and integrating the result on Ij we obtain 


d 

dt 




UjV dx + 



a 


dU^ 

dx 


V dx 


0, W e Vp. 


( 2 ) 


Here, Vp is a finite dimensional space of polynomials of degree up to p. Transforming [xj,Xjj^i] to 
the canonical element [—1,1] by a linear mapping 

xiO = (3) 

yields 

Integrating the flux term by parts we obtain, 

u,vdi + ^ |f/j(i)i''(i) - i7,_.(i)v(-i)] - ^ u,^d( = 0, vv € r,. 


Note that we have chosen an upwind flux so that the value of U at the cell interface is chosen to 
be the value on the left side of the discontinuity. Integrating by parts once more we can write, 

where [[Uj]] = Uj{—1) — denotes the jump between the endpoints of the numerical solution 

at the interface of the j-th and {j — l)-th cells. Note that the second integral in this expression is 
entirely local, as opposed to the term in (|1]). Next, we choose the Legendre polynomials as the basis 
for the finite element space Vp. Recall |T], that the Legendre polynomials T’fc(O) ^ = 0,1, 2,..., 
form an orthogonal system on [—1,1] 


/ VkPidI = — - dki, 

^ 2k +1 ’ 


( 6 ) 


where is the Kroneker delta. With the chosen normalization ([2]), the values of the basis functions 
at the end points of the interval [—1,1] are pQ 

Pt(l) = 1, U(-l) = (-1)‘. (7) 


The numerical solution can be written in terms of this basis as 

p 

Uj = ^ CjkPk, ( 8 ) 

k=0 
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where Cjk is a function of time t. Substituting (|8]) into ([5]), choosing V = Pk, k = 0,1,... ,p, and 
using ([7]) and (jH]) results in the system of equations 


d {2k + l)a 
dt ^ hj 



(2fc+l)a (_i),[|^^ll, 


0 


(9) 


ioT k = 0,1,... ,p. This system of equations is enough to determine Uj dehned by (jH]). However, 
since we are interested in the equation which Uj itself satishes we can reconstruct Uj by multiplying 
each equation in the system by Pk and summing over all k. We can thereby write the following 


exact expression for ^Uj, 





Since the Legendre polynomials are an orthogonal family it is clear that the hrst summed term in 
this expression is simply the projection of into the hnite element space 14- Moreover, since 
is already in the hnite element space this projection is exact. Hence we can write, 

fP + T,Fp = -y 

To simplify this expression further let us use the following proposition: 

Proposition 1. 



where Rpj^i is the right Radau polynomial J7^ of degree p + 1, which is defined as R 
Proof. It is known that the Legendre polynomials satisfy [1], 




Therefore, a simple calculation shows 


- Pv] 

= {-lY[{2p +1)P„+ (2p - 3)Fp_2 + ... - (2p - l)Pp_i - (2p - 5)Fp_3 - ...] 


P 


= 5^(-l)'=(2fc + l)P», 


k=0 

which completes the proof. 


□ 


Using this proposition, we rewrite flTT]) compactly as 



( 12 ) 
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This is an equation which the polynomial approximation Uj will satisfy exactly on the cell Ij. In 
particular, Uj solves the same advection equation as the exact solution u, except with a forcing 
term. Note that when approximating smooth solutions the polynomial approximation Uj will be 
locally order p+1 accurate to the exact solution and thus the jump term at the cell interface [[Uj]] 
will be of order as well. This implies that the solution to flT^ and, hence, the numerical 
approximation will in a sense be close to a solution of the advection equation since the forcing term 
is small. 


3 Fourier Analysis 


To investigate the properties of solutions of flT^ we look at a single Fourier mode solution of the 
form Uj{^, t) = Uj{^, where Uj is a polynomial of degree p in As is customary in Fourier 

analysis, we also assume periodic boundary conditions on /. Using these assumptions on the form 
of Uj{^,t) ill we have that this Fourier mode satisfies the ODE 

+ = (13) 

We can solve this ODE explicitly to obtain that 


U,(e,a;) = U,(-l,n;)eV«+^) + 


.l)P+i 


m]] 


'-1 




ds. 


( 14 ) 


The general solution fITT)) is not necessarily a polynomial in Therefore, fITT)) is too general for 
our purposes. Below, we will look for additional restrictions which ensure that this solution is 
polynomial in We state two lemmas which will help us to rewrite and investigate the integral 
term in fflTl) . 

Lemma 1. The integral term in flTTD satisfies the following relation 


(_l)p+i 




ds 


'p+i' 


u)h,s , _ 

ds = —+ 


{uhfiP+^ 


g{uhj)e 2 


tuns , 

^(C+i) _ 


fi.^hj,i)), ( 15 ) 


where giojhj) is a polynomial of degree p in ujhj and f{u:hj,ff) is a polynomial of degree p in both 
uhj and 

Proof. We begin by integrating the integral in fll4p by parts 




/-I 


i-) 

■ .i)P+i 


continue integrating by parts until the remaining integral vanishes 

/ .- V .V, I 1 




tn obtain 
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where 


f—i '(P+i 

9(^h,) = ( 17 ) 

fe=l 

f — l 'iP+l rjk 

n^hj,Q = ( 18 ) 

Note that g{ujhj) is a polynomial of degree p in uhj, and f{u;hj,^) is polynomial of degree p in uhj 

and Therefore, dehning g{uhj) = + giujhj) in ([TB]) we will obtain eqnation flT^ which 

completes the proof. □ 

Lemma 2. The integral term in flTTD also satisfies 


bThl e={iK-)Tj}-_^^(s) * = _e^K+i) + 


+ 






Proof. We prove this lemma in a similar manner to Lemma 1, i.e., by integrating the integral in 
flTdll by parts, this time in reverse order 




+ 




'-1 


Here, from the dehnition of the Radan polynomial Rp+i, we have nsed 1) = 2(—1)^. We 

continne integrating by parts to obtain 


.l)p+i fi d , (-1)^^^ f ^ 


'-1 




R;f\ ’io 


+ 


(_1)P+1 /^h 


i i I?-.(-2) 


Rp-4.rgo + ---, (20) 


where we dehne Rp+i to be the repeated integrals of the right Radan polynomial, i.e., Rp^i\f) = 
Rp+i(0 and 




= / Rp-i-'^(s) ds. 


p+i 


’-1 


( 21 ) 


Finally, nsing this dehnition with the Canchy integration formnla we can write the polynomials 


Rp+\ as 


RdPiO = 




O’+i {k 

which, when nsed in fl20l) . yields flT^ and completes the proof. 


( 22 ) 

□ 
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From these two lemmas we can establish a useful result regarding the polynomials g and /. 

Corollary 1. 


g{ujhj) 2g{uhY 



( 23 ) 


and, in particular, 


g{ujhj) 


= + C>((a;hj)2p+2), 


(24) 


i.e. is the Fade approximant of . 

Before we state the proof of this corollary let us briefly recall the dehnition of a Bade approximant 

In- 


Definition 1. Given integers m and n and a sufficiently smooth function F{z), the ^ Fade ap¬ 
proximant of F{z) is a rational function where F{z) and Q{z) are polynomials of degree m and 
n, respectively, and satisfy 


P{z) 

QY) 


F(z) + C>(z"*+"+^). 


This Fade approximant is unique up to a constant multiple of the numerator and denominator. It 
is conventional to take Q(0) = 1 so that the Fade approximant is uniquely defined. 


We now proceed to prove the Corollary. 

Froof of Corollary 1. Equating the right hand sides of (ITH]) and (IT^ we obtain 
1 




+ 


{uhjY^^ 


g{uhj)e 


^«+i) _ 


fiujhj,^) =-6 2 




+ 


(_l)P+i 


Rp+liO + ^ 

k=l 


r;F'(o 


Solving this expression for immediately yields fl23|l . Subsequently, evaluating fl2^ at ^ = 1 

we obtain 


/M.,1) ^ _ {-lY^Yooh.Y^^ 

g{ujhj) ^giuhj) 



From the dehnition of Rf^i in terms of the Legendre polynomials we have that i?“^_]^(l) = 0. Fur¬ 
thermore, from the dehnition of (^) in fl22l) and the orthogonality of Rfj^i to all polynomials 

of degree p — 1 we have that = 0 for /c = 1,... ,p. We therefore hnd that 


g{ujhj) 2g{ujhj) 
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which yields 


+ 0{{uh,fP+^). 

g{uhj) 

By Lemma 1, f{z, 1) is a polynomial of degree p while g{z) is a polynomial of degree p + 1 with 
the form g{z) = z^~^^ + g{z). Therefore, we hnd after a possible rescaling that the rational fnnction 
approximates to order 2p+ 2. Therefore it is the uniqne Bade approximant of e^. □ 

Using Lemmas 1 and 2 we can write the general solntion (lT4ll in two ways: 




im] 

iuh,)P+^ 


g{uhj)e - f{ujhj,0 


= Uj-i{l,u)e [9{ujhj)e 2 '- f{uhj,0] , 


(25) 


and 




[Kl] 




k=l 


= £;,_,(i,u.)e^«+‘i +frolic,]) 




k=l 


(0 


P +1 (0 

(26) 


Note that the solntion corresponding to the exact advection of the downwind point Uj-i{l,u) is 

and, hence, from fl25l) and fl26ll the general solution for the numerical approxi¬ 
mation in cell Ij consists of two parts: exact advection of the downwind value in cell Ij-i and higher 
order error terms which are proportional to the magnitude of the jump at that interface. 

As mentioned above, we are interested only in polynomial solutions of flT^ . The reason for this 
is that we know the numerical solution is polynomial in ^ for all times t. Hence, the numerical 
solution should be composed solely of solutions of (IT^ which are polynomials in By examining 
(]25ll we see that the solutions Uj will be polynomial in ^ only when 


U,_i(l,a;) + [[U,]] 


{uh,)P+^ 


0 


(27) 


is satished. Hence, assuming g{uhj) 7 ^ 0, we obtain after rearranging that Uj{—l,u) is related to 
Uj_i{l,uj) by 




gjuhj) - {uhj)P+^ 
g{uhj) 


Using the above relation in fl2^ we obtain after rearranging that the polynomial solutions of flT^ 
have the form 




uj) 


g{uhj) 


(28) 


Thus, we obtain that the polynomial solutions on each cell are completely determined by the rational 
function and the value of the numerical solution at the downwind point of the previous cell. 
















( 29 ) 


Using fl2^ from Corollary 1 we have that at ^ = 1 the polynomial solutions fl28|) satisfy 


(j2^ means that the value of the solution at the right downwind point of cell Ij-i is advected to the 
downwind point in the cell Ij with local accuracy 2p + 2. This is related to the observed strong order 
2p + 1 superconvergence at the downwind point of the numerical solution [2]. In fact, this is also 
related to the study of the superaccurate errors in disipation and dispersion of the DG scheme. The 
same Fade approximant was studied by Hu and Atkins in |9], Ainsworth in [3], and Krivodonova and 
Qin in 110]. In each paper the authors note that the superaccuracies in dissipation and dispersion 
errors stem from the accuracy of this Fade approximant. A key difference here, however, is that we 
have not made the assumption of a uniform mesh. Hence we immediately obtain from fl2^ that we 
can extend the previously known results concerning the local 2p + 2 order of accuracy in dissipation 
and 2p + 3 order of accuracy in dispersion of the DG method to non-uniform meshes. 

In O E] the Fade approximant in was used to relate the numerical frequencies and wavenum¬ 
bers of the scheme by assuming that u was an exact frequency and hnding the numerical wave 
number kn- Another approach was taken in [10] where the authors were interested in the spectrum 
of the DG method, i.e. the precise values of the numerical frequencies u with a hxed wavenumber k. 
In this work we will use the latter approach of estimating the numerical frequencies u for problems 
with periodic boundary conditions. To this end we use the hrst line of fl2^ together with periodicity 
of the boundary conditions to obtain the following condition on u 


TT 1) 

g(uhj) 


(30) 


Hence, the numerical frequency u must satisfy this relation. Solving fl30|) for every value of u, 
however, is difficult since it would require hnding the roots of a high order polynomial. An attempt 
to describe these values for some particular meshes was made in HU, but this is beyond the scope of 
this paper. We will instead make the simplifying assumption that the mesh is uniform and obtain 
that the values of oj are solutions of 


g{ujh) 


(31) 


where is an A^-th root of unity, i.e. Kn = n = 0,..., — 1, where L is the length of 

the domain /. Note that since the mesh is uniform and each downwind point of this solution is 
related by Uj{l,uj) = = e^’*^[/j_i(l,cu) the exact physical frequency for this wave 

is a; = Kn- In the following theorem we give an estimate on the values of u which are solutions of 

(EU). 


Theorem 1. For a uniform mesh of N elements, there are {p + 1)N solutions of flT^ of the 
form Uj{f,t) = Uj{f,u)e~°'‘^^ which are polynomials in f. These solutions satisfy Uj{l,oj) = 
e'^"^Uj_i{l,u) for each j where Kn = n = 0,..., A^ — 1 . For each Kn there are p + 1 values 
u = uo,u)i,... ,LJp which have the expansions 


COo = Kn + 0{KlP+^h^P+^) 
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and 

UJm=^ + 0{Kn), m=l,...,p, 

h 

where prn o^re the p non-zero roots of the polynomial g{z) — f{z) and satisfy Re(/im) > 0. 

Proof. We begin by noting that from 

/M l) 

g{ujh) 

there should be at least one solution of of the form co = Kn + The condition 

flMll itself for the numerical frequency u can be rearranged to obtain 

g{ujh)e^-^ - fiooh, 1) = 0. (32) 

This expression is a polynomial of degree p + 1 in a; and, therefore, has up to p + 1 distinct roots. 

Regarding h as a small parameter, we have from the form of fl32p that we can asymptotically 

approximate each root using the expansion 

d-i , , , 

LO — —-h Uq + Uih + . . . . 

h 

Using this expansion in ([32]), and expanding = 1 + Knh + 0{{Knh)‘^) we obtain 

g{d-i) - /(d-i, 1) + g{d-i)Knh + g'{d-i)doh - /'(d-i, l)doh + 0{{KnhY) = 0. (33) 

Setting the powers of h equal to zero we find that we can determine the leading order asymptotic 
behaviour of each root by finding the possible values of the coefficient d_i which solve 

p(d_i) -/(d_i,l). (34) 

Firstly, evaluating fl 2 T|) at ojhj = 0 gives that /(0,1) = p(0) so d-i = 0 is a root of flMD . Furthermore, 
differentiating ( 12 T|) and evaluating at uhj = 0 yields that /'( 0 , 1 ) ^ g'{0) and, hence, d_i = 0 is 
a simple root. Finally, when this Fade approximant was studied in m, the authors showed that 
non-zero roots of the polynomial g{z) — f{z, 1) lay in the right-half complex plane. Therefore we 
can conclude that there are p roots of the form 

UJm, , (Dy Kn ) , 

h 

where Re(pm) > 0, and one root which corresponds to d-i = 0. Clearly the choice of d-i = 0 must 
correspond to the solution ojq = Kn + thus we obtain the result. □ 

From this Theorem we have that for each there are p-|-1 independent polynomial solutions of 
firsj) which satisfy Uj{l,uj) = for all j. One corresponds to coo = Kn + 

and can be seen as ‘physical’ as it propagates with a numerical frequency which is close to the 
exact frequency. The other, ‘non-physical’, solutions are dampened out exponentially quickly. This 
property of the numerical frequencies of the DG method was conjectured by Guo et al in [S] where 
the authors explicitly calculated similar expansions of the numerical frequencies ojm for F = 1 ; 2 and 
3. Since the non-physical modes are damped out exponentially quickly we see that after sufficiently 
long times the accuracy of the numerical solution will be completely determined by the accuracy 
of the physical mode. Hence, if we specifically choose the initial projection of the exact solution to 
ensure that the physical mode is high-order accurate, we should preserve this accuracy for t > 0 . 
We formalize this observation in the following theorem. 


10 



Theorem 2. Let u{x,t) be a smooth exact solution of ([T]) on the interval I with periodic boundary 
conditions. Let Uh be the numerical solution of the DG scheme on a uniform mesh of N elements 
and let Uj be the restriction of the numerical solution to the cell Ij. Let ej{f,t) = Uj — Uj be the 
numerical error on Ij (mapped to the canonical element [— 1 , l]j. Suppose the projection of the 
initial profile u{x, 0) into the finite element space is chosen such that 

0) - 0)] di = A: = 0,... ,p, (35) 

is satisfied. Then the error on cell Ij will tend exponentially quickly towards the form 

^— [[Uj]]Rp+i + Oip+2{t)Rp+i ^^ + ... + a2p+i{t)Rpj^i (36) 

where ak{t) = 0{h^) and, in particular, 

e,{l,t) = 0{h^^^^). 


Proof. We begin by assuming for simplicity that the exact solution can be written as the sum 

N-l 

u{x,t) = 

n =0 

where Kn = and L is the length of I. The coefficients Un are found by the discrete Fourier 
transform and satisfy 

N-l 

u{xj,0) = j = l,...,N. 

n=0 

Of course, in general the exact solution cannot be written in such a way but provided u is sufficiently 
smooth and N is sufficiently large the error in such an approximation should be negligible compared 
to the error in the polynomial approximation on each cell. Without loss of generality, let us consider 
the numerical approximation of just one of these Fourier modes, u{x,t) = Restricting 

this Fourier mode to the cell Ij and mapping to the canonical element we see that 

Since the mesh is uniform, the projection of «j(a;,0) into the hnite element space will be of the 
form [/j(^,0) = Une'^"^^U{f) for every j and we immediately obtain that Uj{l,t) = e''"^L(,_i(l,t) 
for every j. We therefore can express the numerical solution as a sum of the p + 1 independent 
polynomial solutions found in Theorem 1 which satisfy Uj{f,t) = Lfj{f,ujm)e~°'^""^ and Lfj{l,0Jm) = 
e'^"^17j_i(l, Wm), where the are the p + 1 distinct values which satisfy Hence 


V 

m=0 


f 

g{ujmh) 


Since the physical frequency uq is an accurate approximation of the exact frequency Kn to order 
0{KlP+^h‘^P+^), we have by Corollary 1 the expansion 


^ f{^oh,0 




(-l)P+nwoh)P+^ 

2g{uJoh) 



+0{{KnhYP+^). 
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Therefore performing the initial projection and using this expansion together with the orthogonality 
of the Radau polynomial Rp+i we hnd that 






[U,{^,0)-u,{C,0)]PkdC= I 

_1 g{umh) 


m=l d-1 


Pk dC + 0{KlP+^h^P+^-^). 


Thus, the requirement of the initial projection to satisfy fl3^ will be satished by the choice of Cq = 
Un + 0{h'^P~^^) and X]m=i ^m '^g‘{Yh) ^ where 7 = Hence, this initial projection yields 

a high-order accurate physical mode of the numerical solution. Finally, we know from Theorem 1 
that uji,... ,ujp have positive real parts of order O Hence these components of the solution are 
damped out exponentially quickly in time and the numerical solution tends to the form 

_ ^ KnXj—aUJQt f{^oh,0 




g{uoh) 

We can write this solution in the form (l26l) to hnd that 


+ C>(h 2 p+i). 


(37) 






k=l 


V+l 


(0 


+ C>(h 2 p+ 1 ), 


= %(«.«)+ bd?^iRii 


Rp+iiO + Y 


00 / , \ A: 


k=l 


+ C>(h2p+i) 


From this, we see from that the error for the numerical approximation has the form 

k=l V 2 / 


£:;(«. i) = 2 Jhh|Rll 


Rp+liO + Y 


'P+1 




+ C>(h 2 p+ 1 ), 


and since this is true for any Fourier node we obtain the result by summing this expression 

over all Fourier modes. □ 


From (136|) we see that the numerical solution will be one order more accurate at points ^0 
such that = 0, i.e. the roots of the right Radau polynomial. This is refered to as the 

spatial super convergence of the method. From the proof of Theorem 2 we see that in order for the 
numerical error to tend to a superconvergent form the key requirement of the initial projection is 
that it projects onto the physical mode with high order accuracy. For example, an initial projection 
which consists of simply interpolating the initial data at equidistant points will not satisfy the 
conditions of Theorem 2 and thus we do not observe super convergence of the numerical solution at 
the downwind points at any time. 

Examining the superconvergent form we establish some useful corollaries. First, we note 
that once the non-physical modes have been damped out the remaining physical modes will be 
advected with order accuracy. Hence, for the physical modes, the DG method can be viewed 
as an order 2p -|- 1 scheme. Second, since the initial projection in Theorem 2 produces a high-order 
accurate physical mode, and due to the orthogonality properties of the Radau polynomials, we also 
obtain high-order accuracy of the moments of the numerical solution. We states the results formally 
below. 
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Corollary 2. The accumulation error of the superconvergent numerical solution fl37|) is of order 
2p + 1. That is, after sufficiently long time that the non-physical modes of the numerical solution 
have been damped out, the numerical solution satisfies 

||C/j+i(5.i + oA) - = 

Corollary 3. The superconvergent form of the numerical solution (1361) has the property that the 
m-th moment of the error is order 2p-\-1 — m, i.e. 

In the next section we perform several numerical test to confirm the results of Theorem 2 and 
Corollaries 2 and 3. 


4 Numerical Examples 

In this section we will perform several numerical experiments to conhrm the super convergence 
properties stated in the section above for the DG method for the linear advection equation. Specif¬ 
ically, we will conhrm that on a uniform mesh the numerical solution of the DG method with a 
non-superconvergent initial projection will tend exponentially quickly towards the superconvergent 
form fl36l) . Moreover, we will show that when t is sufficiently large the superconvergent numerical 
solution is advected at order 0{h‘^^~^^). Finally, we will show that the moments of the numerical 
error are also high order accurate after sufficiently long time t. 

Our numerical studies were done on the initial value problem 


u{x^ 0) 
u{-l,t) 

= 0, —1 ^ X < 1, 

= Uo{x), 

= w(lG), 

t > 0, 

(38) 


Uo{x) = sindTTx. 


(39) 


All test below are calculated using an RK-4 time-stepping scheme and a CFL number of to 
minimize the error incurred in time integration. 

Superconvergence from more general initial projections 

In the proof of Theorem 1 we showed that the non-physical waves are damped out like e~~^. 
We therefore expect to observe that a numerical solution with an initial projection satisfying the 
conditions of Theorem 2, to have converged to the superconvergent form (1361) when 

e-T = (9(/i2^'+i), 

where pmm is the non-physical numerical frequency with the smallest real part. Therefore, we 
expect that the numerical solution will be superconvergent when 

t = -‘^P±lo{h\ogh). 
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(a) p = 1 


{h)p = 2 



(c) p = 3 


Figure 1: Semi-log plots of norm of the point-wise error at the downwind points of the numerical 
solution with initial projection as a function of time. Solutions are calculated for the linear 
advection, fl38|l - fl3^ on a uniform mesh of = 64 elements, CFL = 


We can estimate the smallest real part of the non-physical numerical frequencies by explicitly 
calculating the roots of the polynomial g{z) — f{z, 1) and hnding the root with the smallest non¬ 
zero real part. This calculation for p = 1,2, 3, and 4 yields pmin = 6, 3,0.42, and 0.058, respectively. 
Therefore, we see that the smallest real part of the non-physical numerical frequencies is decreasing 
very rapidly as the order p increases. Hence we expect that it will take signihcantly longer for the 
non-physical modes to be damped out as the order of the DG method increases. 

In Figure [U we show the error at the downwind point of the numerical solution as a function of 
time for the p = 1,2, and 3 schemes on a uniform mesh of = 64 elements with the usual initial 
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Projection 

Left Radau Projection 

N 

\m 

r 

ko 

r 

||E|| 

r 

Ik-oll 

r 

16 

7.02e-02 

- 

6 .66e-02 

- 

9.63e-02 

- 

1 .22e-01 

- 

32 

8.40e-03 

3.06 

8.90e-03 

2.90 

1 .22e-02 

2.98 

1 .68e-02 

2.86 

64 

1.04e-03 

3.01 

1.08e-03 

3.04 

1.54e-03 

2.99 

2.13e-03 

2.98 

128 

1.30e-04 

3.00 

1.34e-04 

3.01 

1.93e-04 

2.99 

2.67e-04 

3.00 

256 

1.63e-05 

2.99 

1.67e-05 

3.00 

2.43e-05 

2.99 

3.33e-05 

3.00 


Table 1: Linear advection, fl38|l - fl39|) with p = 1 and with the and left Radau initial projec¬ 
tions. error of the downwind points ||i?|| and of the cell averages ||£o|| are shown together with 
convergence rates, r. Errors are calculated aX t = h. 



Projection 

Left Radau Projection 

N 

\m 

r 

ko 

r 

||E|| 

r 

ko 

r 

16 

5.87e-03 

- 

7.96e-03 

- 

6.65e-03 

- 

7.66e-03 

- 

32 

l.lOe-04 

5.72 

1.86e-04 

5.42 

1.38e-04 

5.59 

2.20e-04 

5.12 

64 

2.74e-06 

5.34 

4.04e-06 

5.52 

3.57e-06 

5.27 

5.54e-06 

5.31 

128 

8.01e-08 

5.10 

l.lOe-07 

5.20 

1.06e-07 

5.07 

1.60e-07 

5.12 

256 

2.47e-09 

5.01 

3.28e-09 

5.07 

3.31e-09 

5.00 

4.87e-09 

5.03 


Table 2: Linear advection, fl38|) - fl39|) with p = 2 and with the and left Radau initial projec¬ 
tions. error of the downwind points ||i?|| and of the cell averages ||£o|| are shown together with 
convergence rates, r. Errors are calculated at f = 4h. 


projection. The error at the downwind point is calculated using the norm of the point-wise 
numerical errors at the downwind points, i.e. ||i?|| = h'^j \Uj{l,t) — Uj{l,t)\. We notice from the 
linear shape of the semi-log plots that the error at the downwind point decays exponentially up to 
some critical time, at which point the error remains relatively constant. We also notice that due to 
the scaling of Pmin h takes signihcantly longer for the error at the downwind points to reach this 
critical time as p increases. In the following numerical test we show that once this critical time is 
reached the error at the downwind points is 

In Tables [T]l3] we show the results of our convergence test for p = 1,2, and 3. In each table we 
present the errors at the downwind points of the cells HEH, and the norm of the numerical 
errors in the cell averages, calculated as 


N 

^oll = 

i=i 

The errors are calculated at f = h,4h, and 35h for the p = 1,2, and 3 methods, respectively, in 
order to allow sufficient time for the non-physical modes to dampen out. We calculate these errors 
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Projection 

Left Radau Projection 

N 


r 

ko 

r 

||E|| 

r 

Ik-oll 

r 

16 

5.14e-04 

- 

1.05e-03 

- 

5.14e-04 

- 

1.06e-03 

- 

32 

2.36e-06 

7.76 

4.39e-06 

7.90 

2.30e-06 

7.80 

4.39e-06 

7.91 

64 

9.17e-09 

8.00 

1.77e-08 

7.95 

9.09e-09 

7.99 

1.82e-08 

7.91 

128 

3.63e-ll 

7.97 

6.93e-ll 

8.00 

3.53e-ll 

8.00 

7.13e-ll 

8.00 

256 

2.75e-13 

7.05 

6.53e-13 

6.73 

2.99e-13 

6.88 

5.83e-13 

6.93 


Table 3: Linear advection, fl38|l - fl39|) with p = 3 and with the and left Radau initial projec¬ 
tions. error of the downwind points ||i?|| and of the cell averages ||£o|| are shown together with 
convergence rates, r. Errors are calculated at f = 35h. 


N 

cT 

1 

o' 

r 

\\U{x,2)-U{xA)\\ 

r 

16 

9.16e-03 

- 

6.59e-03 

- 

32 

2.34e-03 

1.96 

8.34e-04 

2.98 

64 

5.90e-04 

1.99 

1.05e-04 

3.00 

128 

1.48e-04 

2.00 

1.31e-05 

3.00 


Table 4: Linear advection, (I38|) - (l3^ with p = 1 and initial projection. norms of difference in 
numerical solutions at different times. Differences are measured between Uj initially and at f = 2, 
after one period, then between Uj at t = 2 and t = 4, after an additional period. 


for two different initial projections. The hrst is the usual projection while the second is a left 
Radau-like projection which is dehned by 



Uj)Pk 


0 , k 


0,...,p- 1, 


and 

f/,(-l,0)=«j(-l,0). 

These projections, while satisfying the conditions of Theorem 2, are far from the superconvergent 
form fl5BD which can be viewed as close to a right Radau projection of the exact solution. In each 
table we observe the expected 2p -|- 1 rate of convergence in both the error at the downwind points 
of the cells and in the cell averages. 

Order 2p -|- 1 advection of superconvergent solution 

Next, we show that once the non-physical modes of the numerical solution have been damped out, 
the remaining modes are advected at order 2p + 1. To show this we use the initial projection 
and calculate the norm of the difference between numerical solutions after 0,1, and 2 periods. That 
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N 

cT 

1 

o' 

r 

||17(x,2)-f/(x,4)|| 

r 

16 

2.87e-04 

- 

1.03e-05 

- 

32 

3.57e-05 

3.00 

3.24e-07 

4.99 

64 

4.46e-06 

3.00 

l.Ole-08 

5.00 

128 

5.58e-07 

3.00 

3.17e-10 

5.00 


Table 5: Linear advection, (l38ll - (l3^ with p = 2 and initial projection. norms of difference in 
numerical solutions at different times. Differences are measured between Uj initially and at f = 2, 
after one period, then between Uj at t = 2 and t = A, after an additional period. 



Projection 

Left Radau Projection 

N 

dII 

r 

11^211 

r 

IIdII 

r 

11^211 

r 

16 

2.92e-03 

- 

8.27e-03 

- 

3.24e-03 

- 

8.06e-03 

- 

32 

1.12e-04 

4.70 

1.04e-03 

2.99 

1.04e-04 

4.96 

1.04e-03 

2.95 

64 

8.09e-06 

3.79 

1.29e-04 

3.01 

7.97e-06 

3.70 

1.29e-04 

3.01 

128 

5.21e-07 

3.96 

1.61e-05 

3.00 

5.19e-07 

3.94 

1.61e-05 

3.00 

256 

3.28e-08 

3.99 

2.00e-06 

3.00 

3.27e-08 

3.99 

2.01e-06 

3.00 


Table 6: Linear advection, fl38|) - (l3^ with p = 2 and with the and left Radau initial projections. 

norms of the hrst and second moments of the numerical error are shown together with convergence 
rates, r. Errors are calculated at t = Ah. 


is, we calculate these differences as 

N 

\\U{x,0)-U{x,2)\\ = hJ2 

i=i 

In Tables 0] and [5] we see that that the difference between the numerical solution initially and after 
one period converges at the usual p + 1 rate. This is expected since the non-physical modes of 
the solution are present initially, and are However, we also see that that the difference 

between the numerical solution after one and two periods converges with order 2p -|- 1. This shows 
that once the non-physical modes of the solution have been damped out, the remaining physical 
modes are advected at order 2p -|- 1. 

Superconvergence of moments 

Finally, we demonstrate the high order accuracy in the moments of the numerical error for p = 2 
in Table [6l We present the norm of the hrst and second moments of the numerical error in each 
cell. The moments are calculated as 

N 

II h ^ ^ 

i=i 


(U, 


Ui 


Pm di 


'-1 


|f/,(^,0)-f/,(^,2)| di. 


'-1 
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The moments are calculated from the numerical solution using the usual initial projection and 
the left Radau-like projection, as above. From this table we see that the m-th moment of the 
numerical error does indeed achieve the predicted order 2p + 1 — m convergence rate. 


5 Discussion 

We have shown that the numerical solution of the DG method is closely related to the Fade 
approximant of the exponential function e^. Indeed, by hnding the Fourier modes of the PDF (IT^ 
which governs the numerical solution we have shown that the polynomial solutions are completely 
described by the value of the solution at the downwind point of the previous cell and the rational 
function . This rational function has a local expansion in hj in terms of the {p + l)-th 

right Radau polynomial and the anti-derivatives of this polynomial. Furthermore, at the downwind 
point of the cell we have that is the Fade approximant of e^. As discussed in [9] 

and [3], the accuracy of this Fade approximant is what gives rise to the high order accuracies 
in both dissipation and dispersion of the DG scheme known as superaccuracy. The expansion of 
this rational function in terms of the right Radau polynomial and its anti-derivatives is what we 
observe to be the super convergence of the numerical solution at the right Radau points and the 
order 2p -|- 1 superconvergence of the downwind point in each cell. Finally, as studied in [10] and 
shown by equation (l30|) . the spectrum of the DG discretization matrix is directly related to this 
Fade approximant. Therefore, these Fourier modes provide a direct connection between the 
three seemingly disparate properties of superaccuracy, super convergence, and the spectrum of the 
DG method. 

We have shown that for a uniform computational mesh of N elements there exist N polynomial 
solutions which can be viewed as physical components of the numerical wave and pN polynomial 
solutions which are non-physical components. Moreover, these non-physical solutions are damped 
out exponentially quickly in time and, therefore, neglecting time integration errors we can conclude 
that the accuracy of the numerical solution for sufficiently large times is completely determined by 
the accuracy of the initial projection of the exact solution onto the physical modes. Furthermore, 
the DG scheme can be viewed as order 2p -|- 1 accurate on these physical solutions. Using this 
result we proved that for a class of initial projections of the exact solution we expect to obtain a 
numerical solution which is super convergent at both the roots of the right Radau polynomial and 
the downwind points of the cell, after sufficiently long times. In particular, there is a class of initial 
projections which do not initially have order accuracy at the downwind point, but will indeed 
obtain this order of accuracy after sufficient time has elapsed. For these projections the points of 
superconvergence will migrate to the roots of the right Radau polynomial exponentially quickly in 
time. 

We intend to investigate whether this approach would be useful in both identifying the super¬ 
convergence and superaccuracy properties of the DG scheme in non-linear and multidimensional 
problems, and in determining the spectral properties of the DG method in these problems. 
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